### Replication Code for Jim Crow paper with Robinson Data ###
# Uses R version 2.15.2

library(foreign)
# version 0.8-52
library(xtable)
# version 1.7-0

rd <- read.csv("datashare.csv")
st <- read.csv("stateLabels.csv",header=FALSE)


# Add state names to the data set #
colnames(st) <- c("num","name")
st$name  <- as.character(st$name)
stname <- rep("Maine",dim(rd)[1])
profbw <- rep(0,dim(rd)[1])
for(i in 1:dim(rd)[1]){
	stname[i] <- st$name[st$num == rd$state[i]]
	if(rd$race[i]==1){profbw[i] <- rd$total[i+1]/sum(rd$total[i:(i+2)])}
	if(rd$race[i]==2){profbw[i] <- rd$total[i]/sum(rd$total[(i-1):(i+1)])}
	if(rd$race[i]==3){profbw[i] <- rd$total[i-1]/sum(rd$total[(i-2):i])} 
}
rd$stname <- stname
rd$profbw <- profbw

# Comparison of Proportion Black (Jim Crow vs non-Jim Crow)
mean(rd$problk[rd$jcrow==1]) # 20% from text
mean(rd$problk[rd$jcrow==0]) # 1.5% from text

# Comparison of Proportion Foreign-Born White (Jim Crow vs non-Jim Crow)
mean(rd$profbw[rd$jcrow==1]) # 3.4% from text
mean(rd$profbw[rd$jcrow==0]) # 17% from text

# List of Jim Crow States (decreasing in problk) #
rd[rd$jcrow==1 & rd$race==1,c("stname","problk","profbw")][sort(rd$problk[rd$jcrow==1&rd$race==1],decreasing=TRUE,index.return=TRUE)$ix,]

# List of non-Jim Crow States (increasing in problk) #
rd[rd$jcrow==0 & rd$race==1,c("stname","problk","profbw")][sort(rd$problk[rd$jcrow==0&rd$race==1],decreasing=TRUE,index.return=TRUE)$ix,]

# Matched State Subset
m <- subset(rd,stname=="Kansas"  | stname=="Wyoming"| stname=="Indiana"| stname=="Nevada")


# State-level variables
s.illit <- tapply(m$illit,m$state,sum)
s.total <- tapply(m$total,m$state,sum)
s.totalblk <- m$total[seq(3,dim(m)[1],3)]
s.proillit <- s.illit/s.total
s.perblk <- tapply(m$perblk,m$state,mean)
s.problk <- tapply(m$problk,m$state,mean)
s.jcrow <- tapply(m$jcrow,m$state,mean)
s.profbw <- tapply(m$profbw,m$state,mean)
s.jcproblk <- s.problk*s.jcrow
s.jcprofbw <- s.profbw*s.jcrow
s.stname <- unique(m$stname)

# Figure 1

#pdf("logfig.pdf")
par(mfrow=c(1,2))
plot(s.jcrow,log(s.proillit) ,type="n",xlim = c(-.6,1.6),ylab="Log(Illiteracy Rate)",xlab="Jim Crow State",main="Ecological Data", xaxt="n",sub="(a)")
axis(side=1, at=c(0,1), labels=c("No","Yes"))
text(x=s.jcrow,y=I(log(s.proillit)),labels = c("IN (all)","KS (all)","WY (all)","NV (all)" ))
segments(x0=0,y0=log(s.proillit[1])-.01,x1=1,y1=log(s.proillit[2])+.01)
segments(x0=0,y0=log(s.proillit[4])-.01,x1=1,y1=log(s.proillit[3])+.01)

plot(m$jcrow[seq(3,12,3)],log(m$illprop[seq(3,12,3)]/m$illprop[seq(1,10,3)]),type="n",xlim = c(-.6,1.6),ylab="Log(Illiteracy Rate)",xlab="Jim Crow State",main="Individual-Level Data",xaxt="n",sub="(b)")
axis(side=1, at=c(0,1), labels=c("No","Yes"))
text(x=m$jcrow[seq(3,12,3)],y=log(m$illprop[seq(3,12,3)]/m$illprop[seq(1,10,3)]),labels = c("IN (black-white)","KS (black-white)","WY (black-white)","NV (black-white)" )   )
segments(x0=0,y0=log(m$illprop[m$stname=="Indiana"& m$race==3]/m$illprop[m$stname=="Indiana"& m$race==1])+.01,x1=1,y1=log(m$illprop[m$stname=="Kansas"& m$race==3]/m$illprop[m$stname=="Kansas"& m$race==1])-.01)
segments(x0=0,y0=log(m$illprop[m$stname=="Nevada"& m$race==3]/m$illprop[m$stname=="Nevada"& m$race==1])+.01,x1=1,y1=log(m$illprop[m$stname=="Wyoming"& m$race==3]/m$illprop[m$stname=="Wyoming"& m$race==1])-.01)
#dev.off()

# true parameter values #

beta0full <- m$illprop[seq(1,10,3)]
beta1full <- m$illprop[seq(2,11,3)]/m$illprop[seq(1,10,3)]
beta2full <- m$illprop[seq(3,12,3)]/m$illprop[seq(1,10,3)]

## combined data log-likelihood function ##

logl <- function(theta,y,x1,x2,yp,t.wf,t.fbwf,t.bf){
b0 <- theta[1]
b1 <- theta[2]
b2 <- theta[3]
lpf <- t.wf*exp(b0) + t.fbwf*exp(b0 + b1) + t.bf*exp(b0 + b2)
retVal <- sum(log(dpois(x=y,lambda=exp(b0 + b1*x1 +b2*x2)))) + log(dpois(yp,lambda=lpf))
return(retVal)
}



### Sampling Distributions (200 observations) ###

## Equal ##
k <- 200
k.b <- rep(round(.333*k),4)
k.fbw <- rep(round(.333*k),4)
k.w <- rep(round(.333*k),4)
S <- 1000
coefmatInd <- matrix(NA,nc=4,nr=S)
coefmat <- matrix(NA,nc=4,nr=S)
ptm <- proc.time()
set.seed(123)
for(s in 1:S){
illsam.w <- rhyper(nn=rep(1,4),m=m$illit[seq(1,10,3)],n=m$total[seq(1,10,3)]-m$illit[seq(1,10,3)],k=k.w) 
illsam.fbw <- rhyper(nn=rep(1,4),m=m$illit[seq(2,11,3)],n=m$total[seq(2,11,3)]-m$illit[seq(2,11,3)],k=k.fbw)
illsam.b <- rhyper(nn=rep(1,4),m=m$illit[seq(3,12,3)],n=m$total[seq(3,12,3)]-m$illit[seq(3,12,3)],k=k.b)
t.w <- m$total[seq(1,10,3)] - k.w
t.fbw <- m$total[seq(2,11,3)] - k.fbw
t.b <- m$total[seq(3,12,3)] - k.b 
ill.w <- m$illit[seq(1,10,3)] - illsam.w
ill.fbw <- m$illit[seq(2,11,3)] - illsam.fbw
ill.b <- m$illit[seq(3,12,3)] - illsam.b
initB0 <- log(tapply(m$illit,m$state,sum)/tapply(m$total,m$state,sum))
for(p in 1:4){
	x1 <- rep(0,k)
	if(k.fbw[p] > 0){x1[(k.w[p]+1):(k.w[p]+k.fbw[p])] <- 1}
	x2 <- rep(0,k)
	if(k.b[p] > 0){ x2[(k.w[p]+k.fbw[p]+1):k] <- 1}
	y <- rep(0,k)
	if(illsam.w[p] > 0){y[1:illsam.w[p]] <- 1}
	if(illsam.fbw[p] > 0){y[(k.w[p]+1):(k.w[p]+illsam.fbw[p])] <- 1}
	if(illsam.b[p] > 0){y[(k.w[p]+k.fbw[p]+1):(k.w[p]+k.fbw[p]+illsam.b[p])] <- 1}
	yp <- sum(c(ill.w[p],ill.fbw[p],ill.b[p]))
	coefmatInd[s,p] <- glm(y~x2+x1,family=poisson)$coef[2]
	badStart = TRUE
	sinit <-0
	while(badStart){
		initPar <- c(log(beta0full[p]), log(beta1full[p]), log(beta2full[p])) + rnorm(n=3,mean=0,sd=sinit)
		if(logl(initPar,y=y,x1=x1,x2=x2,yp=yp,t.wf=t.w[p],t.fbwf=t.fbw[p],t.bf=t.b[p]) > -99999){badStart = FALSE}
	sinit <- sinit + .01	
	}
		topt <- optim(par=initPar, logl, method="BFGS", control=list(fnscale=-1,maxit=500),y=y,x1=x1,x2=x2,yp=yp,t.wf=t.w[p],t.fbwf=t.fbw[p],t.bf=t.b[p])
	coefmat[s,p] <- topt$par[3]
	}
}
proc.time() - ptm

coefmat200eq <- coefmat
coefmatInd200eq <- coefmatInd

## White/Black ##

k <- 200
k.b <- rep(round(.5*k),4)
k.fbw <- rep(0,4)
k.w <- rep(round(.5*k),4)
S <- 1000
coefmatInd <- matrix(NA,nc=4,nr=S)
coefmat <- matrix(NA,nc=4,nr=S)
ptm <- proc.time()
set.seed(123)
for(s in 1:S){
illsam.w <- rhyper(nn=rep(1,4),m=m$illit[seq(1,10,3)],n=m$total[seq(1,10,3)]-m$illit[seq(1,10,3)],k=k.w) 
illsam.fbw <- rhyper(nn=rep(1,4),m=m$illit[seq(2,11,3)],n=m$total[seq(2,11,3)]-m$illit[seq(2,11,3)],k=k.fbw)
illsam.b <- rhyper(nn=rep(1,4),m=m$illit[seq(3,12,3)],n=m$total[seq(3,12,3)]-m$illit[seq(3,12,3)],k=k.b)
t.w <- m$total[seq(1,10,3)] - k.w
t.fbw <- m$total[seq(2,11,3)] - k.fbw
t.b <- m$total[seq(3,12,3)] - k.b 
ill.w <- m$illit[seq(1,10,3)] - illsam.w
ill.fbw <- m$illit[seq(2,11,3)] - illsam.fbw
ill.b <- m$illit[seq(3,12,3)] - illsam.b
initB0 <- log(tapply(m$illit,m$state,sum)/tapply(m$total,m$state,sum))
for(p in 1:4){
	x1 <- rep(0,k)
	if(k.fbw[p] > 0){x1[(k.w[p]+1):(k.w[p]+k.fbw[p])] <- 1}
	x2 <- rep(0,k)
	if(k.b[p] > 0){ x2[(k.w[p]+k.fbw[p]+1):k] <- 1}
	y <- rep(0,k)
	if(illsam.w[p] > 0){y[1:illsam.w[p]] <- 1}
	if(illsam.fbw[p] > 0){y[(k.w[p]+1):(k.w[p]+illsam.fbw[p])] <- 1}
	if(illsam.b[p] > 0){y[(k.w[p]+k.fbw[p]+1):(k.w[p]+k.fbw[p]+illsam.b[p])] <- 1}
	yp <- sum(c(ill.w[p],ill.fbw[p],ill.b[p]))
	coefmatInd[s,p] <- glm(y~x2,family=poisson)$coef[2]
	badStart = TRUE
	sinit <-0
	while(badStart){
		initPar <- c(log(beta0full[p]), log(beta1full[p]), log(beta2full[p])) + rnorm(n=3,mean=0,sd=sinit)
		if(logl(initPar,y=y,x1=x1,x2=x2,yp=yp,t.wf=t.w[p],t.fbwf=t.fbw[p],t.bf=t.b[p]) > -99999){badStart = FALSE}
	sinit <- sinit + .01	
	}
		topt <- optim(par=initPar, logl, method="BFGS", control=list(fnscale=-1,maxit=500),y=y,x1=x1,x2=x2,yp=yp,t.wf=t.w[p],t.fbwf=t.fbw[p],t.bf=t.b[p])
	coefmat[s,p] <- topt$par[3]
	}
}
proc.time() - ptm

coefmat200wb <- coefmat
coefmatInd200wb <- coefmatInd

### Sampling Distributions (300 observations) ###

## Equal ##
k <- 300
k.b <- rep(round(.333*k),4)
k.fbw <- rep(round(.333*k),4)
k.w <- rep(round(.333*k),4)
S <- 1000
coefmatInd <- matrix(NA,nc=4,nr=S)
coefmat <- matrix(NA,nc=4,nr=S)
ptm <- proc.time()
set.seed(123)
for(s in 1:S){
illsam.w <- rhyper(nn=rep(1,4),m=m$illit[seq(1,10,3)],n=m$total[seq(1,10,3)]-m$illit[seq(1,10,3)],k=k.w) 
illsam.fbw <- rhyper(nn=rep(1,4),m=m$illit[seq(2,11,3)],n=m$total[seq(2,11,3)]-m$illit[seq(2,11,3)],k=k.fbw)
illsam.b <- rhyper(nn=rep(1,4),m=m$illit[seq(3,12,3)],n=m$total[seq(3,12,3)]-m$illit[seq(3,12,3)],k=k.b)
t.w <- m$total[seq(1,10,3)] - k.w
t.fbw <- m$total[seq(2,11,3)] - k.fbw
t.b <- m$total[seq(3,12,3)] - k.b 
ill.w <- m$illit[seq(1,10,3)] - illsam.w
ill.fbw <- m$illit[seq(2,11,3)] - illsam.fbw
ill.b <- m$illit[seq(3,12,3)] - illsam.b
initB0 <- log(tapply(m$illit,m$state,sum)/tapply(m$total,m$state,sum))
for(p in 1:4){
	x1 <- rep(0,k)
	if(k.fbw[p] > 0){x1[(k.w[p]+1):(k.w[p]+k.fbw[p])] <- 1}
	x2 <- rep(0,k)
	if(k.b[p] > 0){ x2[(k.w[p]+k.fbw[p]+1):k] <- 1}
	y <- rep(0,k)
	if(illsam.w[p] > 0){y[1:illsam.w[p]] <- 1}
	if(illsam.fbw[p] > 0){y[(k.w[p]+1):(k.w[p]+illsam.fbw[p])] <- 1}
	if(illsam.b[p] > 0){y[(k.w[p]+k.fbw[p]+1):(k.w[p]+k.fbw[p]+illsam.b[p])] <- 1}
	yp <- sum(c(ill.w[p],ill.fbw[p],ill.b[p]))
	coefmatInd[s,p] <- glm(y~x2+x1,family=poisson)$coef[2]
	badStart = TRUE
	sinit <-0
	while(badStart){
		initPar <- c(log(beta0full[p]), log(beta1full[p]), log(beta2full[p])) + rnorm(n=3,mean=0,sd=sinit)
		if(logl(initPar,y=y,x1=x1,x2=x2,yp=yp,t.wf=t.w[p],t.fbwf=t.fbw[p],t.bf=t.b[p]) > -99999){badStart = FALSE}
	sinit <- sinit + .01	
	}
		topt <- optim(par=initPar, logl, method="BFGS", control=list(fnscale=-1,maxit=500),y=y,x1=x1,x2=x2,yp=yp,t.wf=t.w[p],t.fbwf=t.fbw[p],t.bf=t.b[p])
	coefmat[s,p] <- topt$par[3]
	}
}
proc.time() - ptm
#Note: there will be a warning with respect to the fitting of the Poisson model using only the individual level data. This is addressed in the supplementary material.

coefmat300eq <- coefmat
coefmatInd300eq <- coefmatInd

## White/Black ##

k <- 300
k.b <- rep(round(.5*k),4)
k.fbw <- rep(0,4)
k.w <- rep(round(.5*k),4)
S <- 1000
coefmatInd <- matrix(NA,nc=4,nr=S)
coefmat <- matrix(NA,nc=4,nr=S)
ptm <- proc.time()
set.seed(123)
for(s in 1:S){
illsam.w <- rhyper(nn=rep(1,4),m=m$illit[seq(1,10,3)],n=m$total[seq(1,10,3)]-m$illit[seq(1,10,3)],k=k.w) 
illsam.fbw <- rhyper(nn=rep(1,4),m=m$illit[seq(2,11,3)],n=m$total[seq(2,11,3)]-m$illit[seq(2,11,3)],k=k.fbw)
illsam.b <- rhyper(nn=rep(1,4),m=m$illit[seq(3,12,3)],n=m$total[seq(3,12,3)]-m$illit[seq(3,12,3)],k=k.b)
t.w <- m$total[seq(1,10,3)] - k.w
t.fbw <- m$total[seq(2,11,3)] - k.fbw
t.b <- m$total[seq(3,12,3)] - k.b 
ill.w <- m$illit[seq(1,10,3)] - illsam.w
ill.fbw <- m$illit[seq(2,11,3)] - illsam.fbw
ill.b <- m$illit[seq(3,12,3)] - illsam.b
initB0 <- log(tapply(m$illit,m$state,sum)/tapply(m$total,m$state,sum))
for(p in 1:4){
	x1 <- rep(0,k)
	if(k.fbw[p] > 0){x1[(k.w[p]+1):(k.w[p]+k.fbw[p])] <- 1}
	x2 <- rep(0,k)
	if(k.b[p] > 0){ x2[(k.w[p]+k.fbw[p]+1):k] <- 1}
	y <- rep(0,k)
	if(illsam.w[p] > 0){y[1:illsam.w[p]] <- 1}
	if(illsam.fbw[p] > 0){y[(k.w[p]+1):(k.w[p]+illsam.fbw[p])] <- 1}
	if(illsam.b[p] > 0){y[(k.w[p]+k.fbw[p]+1):(k.w[p]+k.fbw[p]+illsam.b[p])] <- 1}
	yp <- sum(c(ill.w[p],ill.fbw[p],ill.b[p]))
	coefmatInd[s,p] <- glm(y~x2,family=poisson)$coef[2]
	badStart = TRUE
	sinit <-0
	while(badStart){
		initPar <- c(log(beta0full[p]), log(beta1full[p]), log(beta2full[p])) + rnorm(n=3,mean=0,sd=sinit)
		if(logl(initPar,y=y,x1=x1,x2=x2,yp=yp,t.wf=t.w[p],t.fbwf=t.fbw[p],t.bf=t.b[p]) > -99999){badStart = FALSE}
	sinit <- sinit + .01	
	}
		topt <- optim(par=initPar, logl, method="BFGS", control=list(fnscale=-1,maxit=500),y=y,x1=x1,x2=x2,yp=yp,t.wf=t.w[p],t.fbwf=t.fbw[p],t.bf=t.b[p])
	coefmat[s,p] <- topt$par[3]
	}
}
proc.time() - ptm
#Note: there will be some warnings with respect to the fitting of the Poisson model using only the individual level data. This is addressed in the supplementary material.  

coefmat300wb <- coefmat
coefmatInd300wb <- coefmatInd


### Table 1 ###

tmat <- rbind(
c(apply(coefmat200eq,2,sd)/apply(coefmatInd200eq,2,sd),apply(coefmat300eq,2,sd)/apply(coefmatInd300eq,2,sd)),
c(apply(coefmat200wb,2,sd)/apply(coefmatInd200wb,2,sd),apply(coefmat300wb,2,sd)/apply(coefmatInd300wb,2,sd)))
xtable(tmat,digits=2)


## Optimal Sampling for Combined Inference ##

# State Numbers
# Indiana 11
# Kansas 21
# Wyoming 41
# Nevada 46

## Conditional Expected Information Function ##
infoFunc <- function(betaw,betac,n00,n10,n01,n11,k){
	n <- n00+n10+n01
	n.nots = n - k
	xbar = (n11 + n10)/n
	zbar = (n11 + n01)/n
	k.xbar.seq <- 0:k
	k.zbar.seq <- 0:k 
	infoMat <- matrix(0,nr=k+1,nc=k+1)
	k11max<- matrix(0,nr=k+1,nc=k+1)
	for(i in 1:(k+1)){ 
		for(j in 1:((k+1)-(i-1))){
			k1. <- k.xbar.seq[i]
			k.1 <- k.zbar.seq[j]
			k11 <- 0
			k10= k.xbar.seq[i];n10.nots = n10-k10
			k01= k.zbar.seq[j] - k11;n01.nots = n01-k01
			k00 = k - (k10+k01) ;n00.nots = n00-k00
			pi.p <- n00 + n10*exp(betaw)+n01*exp(betac)
			pi.p.nots <- n00.nots + n10.nots*exp(betaw)+n01.nots*exp(betac)
			a <- n10.nots*exp(betaw) -(1/pi.p.nots)*(n10.nots*exp(betaw))^2
			b <- n01.nots*exp(betac) -(1/pi.p.nots)*(n01.nots*exp(betac))^2
			d <-  - (1/pi.p.nots)*(n10.nots*exp(betaw))*(n01.nots*exp(betac))
			e <- n10*exp(betaw) -(1/pi.p)*(n10*exp(betaw))^2
			f <- n01*exp(betac) -(1/pi.p)*(n01*exp(betac))^2
			g <-  - (1/pi.p)*(n10*exp(betaw))*(n01*exp(betac))
			Info <- matrix(c(e,g,g,f),nr=2) - matrix(c(a,d,d,b),nr=2)
			infoMat[i,j] <- Info[1,1] - Info[1,2]*Info[2,1]/Info[2,2]
			}
		}	
	return(infoMat)		
	}
##
	
## Numbers for Table 2 ##

# Rows 1 and 2
betaw = 2.5
betac = 2.5
s.seq <- unique(m$state)
state.max <- matrix(0,nr=4,nc=2)
for(s in 1:4){
	infoMat <- infoFunc(betaw=betaw,betac=betac,n00=m$total[m$state==s.seq[s]][1],n10=m$total[m$state==s.seq[s]][3],n01=m$total[m$state==s.seq[s]][2],n11=0,k=100)
	state.max[s,] <- which(infoMat==max(infoMat),arr.ind=TRUE)
	}
sm <-state.max - 1
xtable(t(sm),digits=0)

# Rows 3 and 4
betaw = 0.0
betac = 0.0
s.seq <- unique(m$state)
state.max <- matrix(0,nr=4,nc=2)
for(s in 1:4){
	infoMat <- infoFunc(betaw=betaw,betac=betac,n00=m$total[m$state==s.seq[s]][1],n10=m$total[m$state==s.seq[s]][3],n01=m$total[m$state==s.seq[s]][2],n11=0,k=100)
	state.max[s,] <- which(infoMat==max(infoMat),arr.ind=TRUE)
	}
sm <-state.max - 1
xtable(t(sm),digits=0)

## Additional Sampling Distribution for Table 3 ##

# no restriction #
k <- 300
k.b <- c(floor(.69*k), floor(.71*k), floor(.37*k), floor(.32*k))
k.fbw <- k - k.b  
k.w <- rep(round(0*k),4) 
S <- 1000
coefmat <- matrix(NA,nc=4,nr=S)
ptm <- proc.time()
set.seed(123)
for(s in 1:S){
illsam.w <- rhyper(nn=rep(1,4),m=m$illit[seq(1,10,3)],n=m$total[seq(1,10,3)]-m$illit[seq(1,10,3)],k=k.w) 
illsam.fbw <- rhyper(nn=rep(1,4),m=m$illit[seq(2,11,3)],n=m$total[seq(2,11,3)]-m$illit[seq(2,11,3)],k=k.fbw)
illsam.b <- rhyper(nn=rep(1,4),m=m$illit[seq(3,12,3)],n=m$total[seq(3,12,3)]-m$illit[seq(3,12,3)],k=k.b)
t.w <- m$total[seq(1,10,3)] - k.w
t.fbw <- m$total[seq(2,11,3)] - k.fbw
t.b <- m$total[seq(3,12,3)] - k.b 
ill.w <- m$illit[seq(1,10,3)] - illsam.w
ill.fbw <- m$illit[seq(2,11,3)] - illsam.fbw
ill.b <- m$illit[seq(3,12,3)] - illsam.b
initB0 <- log(tapply(m$illit,m$state,sum)/tapply(m$total,m$state,sum))
for(p in 1:4){
	x1 <- rep(0,k)
	if(k.fbw[p] > 0){x1[(k.w[p]+1):(k.w[p]+k.fbw[p])] <- 1}
	x2 <- rep(0,k)
	if(k.b[p] > 0){ x2[(k.w[p]+k.fbw[p]+1):k] <- 1}
	y <- rep(0,k)
	if(illsam.w[p] > 0){y[1:illsam.w[p]] <- 1}
	if(illsam.fbw[p] > 0){y[(k.w[p]+1):(k.w[p]+illsam.fbw[p])] <- 1}
	if(illsam.b[p] > 0){y[(k.w[p]+k.fbw[p]+1):(k.w[p]+k.fbw[p]+illsam.b[p])] <- 1}
	yp <- sum(c(ill.w[p],ill.fbw[p],ill.b[p]))
	badStart = TRUE
	sinit <-0
	while(badStart){
		initPar <- c(log(beta0full[p]), log(beta1full[p]), log(beta2full[p])) + rnorm(n=3,mean=0,sd=sinit)
		if(logl(initPar,y=y,x1=x1,x2=x2,yp=yp,t.wf=t.w[p],t.fbwf=t.fbw[p],t.bf=t.b[p]) > -99999){badStart = FALSE}
	sinit <- sinit + .01	
	}
		topt <- optim(par=initPar, logl, method="BFGS", control=list(fnscale=-1,maxit=500),y=y,x1=x1,x2=x2,yp=yp,t.wf=t.w[p],t.fbwf=t.fbw[p],t.bf=t.b[p])
	coefmat[s,p] <- topt$par[3]
	}
}
proc.time() - ptm

coefmat300opt <- coefmat

## Table 3 ##
tmat <- rbind(
apply(coefmat300wb,2,sd),
apply(coefmat300opt,2,sd)
)
xtable(tmat,digits=2)

## Table A1 ##
# Indiana #
rbind(
c(m$illit[1],m$total[1],m$illprop[1]*100),
c(m$illit[2],m$total[2],m$illprop[2]*100),
c(m$illit[3],m$total[3],m$illprop[3]*100),
c(sum(m$illit[1:3]),sum(m$total[1:3]),sum(m$illit[1:3])/sum(m$total[1:3])*100 )
)

# Nevada #
rbind(
c(m$illit[10],m$total[10],m$illprop[10]*100),
c(m$illit[11],m$total[11],m$illprop[11]*100),
c(m$illit[12],m$total[12],m$illprop[12]*100),
c(sum(m$illit[10:12]),sum(m$total[10:12]),sum(m$illit[10:12])/sum(m$total[10:12])*100 )
)

## Table A2 ##

# Kansas #
rbind(
c(m$illit[4],m$total[4],m$illprop[4]*100),
c(m$illit[5],m$total[5],m$illprop[5]*100),
c(m$illit[6],m$total[6],m$illprop[6]*100),
c(sum(m$illit[4:6]),sum(m$total[4:6]),sum(m$illit[4:6])/sum(m$total[4:6])*100 )
)

# Wyoming #
rbind(
c(m$illit[7],m$total[7],m$illprop[7]*100),
c(m$illit[8],m$total[8],m$illprop[8]*100),
c(m$illit[9],m$total[9],m$illprop[9]*100),
c(sum(m$illit[7:9]),sum(m$total[7:9]),sum(m$illit[7:9])/sum(m$total[7:9])*100 )
)


